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ABSTRACT 

An elastic large displacement thick-shell mixed finite 
element is modified to allow for the calculation of viscoelastic 
stresses. Internal strain variables are introduced at the 
element's stress nodes and are employed to construct a viscous 
material model. First order ordinary differential equations 
relate the internal strain variables to the corresponding elastic 
strains at the stress nodes. The viscous stresses are computed 
from the internal strain variables using viscous moduli which are 
a fraction of the elastic moduli. The energy dissipated by the 
action of the viscous stresses is included in the mixed 
variational functional. The nonlinear quasi-static viscous 

equilibrium equations are then obtained. Previously developed 

Taylor expansions of the nonlinear elastic equilibrium equations 
are modified to include the viscous terms. A predictor— corrector 
time marching solution algorithm is employed to solve the 
algebraic-differential equations. The viscous shell element is 
employed to computationally simulate a stair-step loading and 
unloading of an aircraft tire in contact with a frictionless 
surface. 
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INTRODUCTION 


tires are composite structures manufactured with 
viscoelastic materials such as carbon black filled rubber and 
nylon cords. When loaded, tires experience large deflections and 
moderately large strains . 1 Finite element models of tires 
typically employ either two-dimensional thick shell or three- 
dimensional solid elements . 2,3 Elastic finite element shell 
models for tires have been used to predict the shape of tire 
footprints as a function of loading . 4-9 Elastic models do not 

include the viscoelastic nature of the tire which can have a 
significant effect on load-displacement curves. This study is the 
computational part of an effort in which the quasi-static 
viscoelastic loading and unloading of an aircraft tire was 
measured and computationally simulated. The experimental effort 
and the load measurements have been reported elsewhere . 10 

In several previous studies viscoelastic constitutive models 
have been utilized to determine the dynamic deformations of tires. 
The following references are provided as a starting point for 
readers interested in obtaining details about other viscoelastic 
finite element models for tires. Padovan, et al . 11,12,13 

performed an extensive study in which a finite element algorithm 
was developed for rolling tires. Padovan 's model included the 
effects of large deformations and contact. It also employed 
fractional derivatives to model the viscoelastic effects. The 

tire models made by Oden et al . 2,3 also included the effects of 

large deformations and contact. However, Oden's model employed 
the history integral formulation for the viscoelastic effects . 

In this report, internal strain variables are employed to 
convert an elastic mixed shell element into a viscous mixed shell 

element. The model is developed as follows. Internal strain 
variables are introduced at the stress nodes of the mixed element. 
First order differential equations relate the internal variables 
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to the physical strain variables . The equations represent a 
Maxwell solid. ' Viscous stresses are determined from the 

internal strains by using material parameters referred to as 
viscous moduli. An expression for the energy dissipated during 
deformation is computed from the viscous stresses. This is 
accomplished by employing the finite element interpolations that 
are used to compute the elastic stresses from the elastic strains 
in the elastic version of the shell element. The dissipation 
energy functional is added to the mixed variational statement for 
the elastic problem. Nonlinear algebraic equilibrium equations 
are determined and are numerically solved, simultaneously, with 
the internal variable differential equations. The numerical 
solution procedure employs the Newton-Raphson method for the 
nonlinear algebraic equations , and the trapezoidal method for the 
differential equations in a predictor-corrector combination. The 

tangent matrix required in the Newton-Raphson scheme is a modified 

4 7-9 

version of the previously determined ' tangent matrix for the 
nonlinear elastic problem. 

At the end of the paper, the viscous shell element is 
employed in a computational simulation of a stair-step loading and 
unloading of an aircraft tire. The elastic material constants 

7 

used are from a previous effort. The viscous material constants 

required are estimated by making a least-squares fit of load- 
relaxation data, as predicted by a one-dimensional version of the 
viscous model, to measured load- relaxation data. In the 

simulation, the stair-step tire rim displacement, as measured in 

the experimental effort, ^ is enforced. The computed and measured 

stair-step hysteresis curves are presented together for reference. 


VISCOELASTIC MIXED SHELL ELEMENT 

An elastic shell element capable of modeling geometrically 
nonlinear deformations of thick laminated composites was developed 
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by Noor, et al. Figure 1. shows the physical variables 

employed to describe the energy in the deformed shell. The 
elastic finite element has nine displacement nodes with five 
variables at each node, and four stress nodes with eight variables 
at each node, see Figure 2. The constitutive model 7-9 in the 
shell generalized coordinates is abbreviated as 


° ~ C e (1) 

where <J — {N s ,Nq,N s q, M s ,Mq,M s q,Q s ,Qq) is a vector of stress 

variables, e = (e s ,£ d ,2e se ,K s ,K e ,2K sd ,2£ s3 ,2e e3 ) T is a vector of the 

19 20 

Sanders -Budiansky ' nonlinear strains, and C is a matrix of 

elastic stiffness constants. The elastic element employs the 
Hsllinger— Resiner mixed variational principle which is constructed 
as follows. The complementary form of the energy is integrated 
over the volume of the shell and the total work done by external 
forces is subtracted. This results in which is expressed as 

follows . 


It hr ~ |(^ Te — W (2) 

where Q is the volume of the shell, F is a flexibility matrix, 
and W is the work done by external forces. 

Next, the energy functional, 71^, discretized by the 

. . 4 

f ini te element method. At the element level, the displacements 

and stress resultants are approximated by employing interpolation 
functions with the nodal values shown in Figure 2 . The Sanders— 
Budiansky nonlinear strains are computed and substituted into 
Equation (2) above. After the volume integration is performed for 
an element, the Hellinger-Resiner variational expression is given 
in short-hand notation as follows. 
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n%{xM) = v-u c -w 


(3) 


where V= \o\dCl = h T (s fa +±M^J = h T e, 

^e// 

u c = f -c T F<Jdn = — h T F h, W= x T p = p T x, x and e are 

J 2 2 

n e , t 

vectors containing the element ' s nodal displacements and strains , 

A A A 

h is an element level vector, and are operators that 

produce the linear and nonlinear contributions to e from the nodal 

* 

displacements, r is an element level flexibility matrix, and p 
represents the consistent applied load vector. 

Internal strain variables are employed herein to modify the 

above formulation making it applicable to a Maxwell type 

14 15 

viscoelastic material. ' Introduce, in vector form, internal 
strain fields, ey v , and matrices of viscous stiffness constants, 
C v , within the element . The viscous stress vector at a point in 
the element is computed as follows. 


<7v=X C v e ,v (4) 

7=1 

The total value of the stress vector, elastic plus viscous, at a 
point is <J t =<J + <J V . Next, following the computation of the 

elastic potential, V, in Equation (3), the energy dissipated by 
the viscous stresses, Q, throughout the element is computed as 
follows . 

P r-p A i-p A rp / A I A \ 

Q= \(J v edCl = h y e = h v ^S & +-M Bte J (5) 

Q-elt 
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The mixed variational 


where h v is an element level vector, 
functional for the viscous element, 7l^ v , is given by 

Khr{**)=V + Q-U c -W (6) 

The model is completed by relating the rate of change of the 
internal strain variables to the physical strain variables. Here 
we employ a simple form of the Maxwell solid theory by reguiring 
the following differential equations to be satisfied. 14 ' 15,18 


( _ de 

dt r jv dt 


j = »,n 


(7) 


Equations (4) and (7) determine the viscous stresses for a time 
dependent deformation of the element, x(f) . Note, an advantage of 

this algorithm is that a variety of viscoelastic models can be 
employed by simply changing Equation ( 7 ) . 

At each instant of time, the element equilibrium equations 

are given by the first variation of 7E ^j v , Equation (6). The 
equilibrium equations are 


f fi (x,h) = (^S & +lM„ ta j-Fh= 0 

and 

f i (x,h) = (s^+M n ^) T (h + h v )-p = 0 


( 8 ) 

(9) 


where and are the derivatives of the operators and 

A 

with respect to the element displacement variables , x . 
Equations (8) and (9) are assembled by standard methods to obtain 
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the global equilibrium equations . The global equations are then 
solved simultaneously with Equation 7 (for all elements.) 

The Taylor expansion of Equations (8) and (9) produces the 
element's tangent matrix. The resulting element level Newton- 

A ^ 

Raphson equations for the increments of the variables Ah and Ax 
are given below. 

-F S, + M n(x 

(S (+ M„ ( ,) r M„,(h + h v ) 



( 10 ) 


where is the second derivative of the operator with 
respect to the nodal displacements. Equations (10) are assembled 
for all elements and solved to provide estimates of the variables 
increments across a time step. The new elastic strains at the end 
of the time step are computed at all stress nodes. The internal 
strain variables are then estimated at the end of the time step by 
employing the trapezoidal method to Equation (7). Next, global 
equilibrium is checked. If equilibrium is not satisfied the 
process is repeated. When equilibrium is satisfied the required 
output is computed and the time is advanced. 


TIRE STAIR-STEP LOADING SIMULATION 

The aircraft tire simulated below is a 32 x 8.8, type VII, 

bias-ply Shuttle nose-gear tire which has a 20-ply rated carcass 

7-10 

and a maximum speed rating of 217 knots. The tread pattern 

consists of three circumferential grooves and the rated inflation 

10 

pressure is 320 psi. During the stair-step test, the tire was 

inflated to 300 psi. The rated operating load for the tire is 

15,000 pounds. 

Elastic Material Model 
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are 


The details of the tire's elastic material model 
described by Tanner . The tire is a cord-rubber composite and was 

treated as a laminated material. It was divided into seven 
regions in the direction of the meridian (from the center of the 
tread region to the rim. ) Tire thickness, properties of the 
plies, etc. were measured and tabulated. Elastic constants were 
computed by the law of mixtures to obtain linear orthotropic 
stress-strain constitutive models for each layer. These 

properties were transformed to the shell coordinate system and 
integrated through the thickness of the shell elements. 

Viscoela stic Material Model 

One dimensional tire loading data 10 is employed below to 

°ktain an approximation to the viscous material properties . 
Obtaining and analyzing time dependent multi-axial test data would 
result in a more accurate viscous material model than the model 
described in this study. The task of performing multi-axial 
stress— strain tests on coupons cut from the tire was beyond the 
scope of this effort. Below we describe how the measured one- 
dimensional load-relaxation data is least— squares fit to a 
simplified form of Equation (4). 

When a material is rapidly deformed from a relaxed unstrained 
state into a deformed state in which the strains are held 
constant, we have a state of stress relaxation in the material. 
In this case. Equations (7) can be integrated. The solutions are 
substituted into Equation (4) and the following stress-relaxation 
equation for G t (t ) is obtained. 




N 


-t 


C+C,x«, c '' 

rt=l 



( 11 ) 


where e(? — 0 ) is the initial value of the rapidly obtained strain 
state achieved. 


8 



Since the test data is one dimensional we can not determine 
the matrix C v . Here, C v is assumed to be equal to the matrix C 
throughout the tire. With this assumption, Equation (11) becomes 




N V 

1+ Xa„e • 


n = 1 


C e(f = 0 + ) 


( 12 ) 


Note, modeling the tire with shell elements neglects three- 
dimensional effects (such as tread compression, bead-rim 
interaction, etc.) When this shell element tire model is suddenly 
pressed against a frictionless platform. Equation (12) implies 
that the total measured load on the platform will be given by 



l + £a n e T ' 

n—\ 


f(‘ = 0 + ) 


(13) 


where /(( = 0 + ) is the elastic part of the tire's load on the 

platform. The Prony series coefficients in Equation (13) are 
determined below by performing a constrained least-squares fit to 
measured tire load-relaxation data. 

The classical procedure of inspecting the log (load) versus 
log (time) data curve for a load-relaxation test is used to select 

a range of time constants, {t„} The range selected must 

include the full spectrum of decay rates needed to model the 
relaxation curve. The error, at time t { is 




n y- 

1+ I a n e n 

n = 1 


/(' = 0+ ) 


(14) 
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where /(*,•) is the measured platform load at time t.. The square of 

the error indicated by Equation (14) is added for all the data and 
is minimized with the constraint that the coefficients, a n , be 

positive. That is, so that > o}^ =1 . The algorithm employed to 

solve this constrained minimization problem is due to Rusin. It 

22 

was also used by Johnson and Quigley as part of an algorithm 

which solves nonlinear frictionless contact problems. 

The relaxation data shown in Figure 3 is from a 900 sec load- 
relaxation test. In the test, a large displacement of the tire's 
rim was rapidly enforced and then held approximately constant for 

900 sec. The enforced displacement was selected 10 so that the 

total tire loading would be near 20,000 lbs. Inspection of the 
tire rim to platform displacement data indicated that the platform 
was moving slightly during the test. The elastic component of the 
load-relaxation data was adjusted, as described below, to correct 
for the platform drift. 

The adjustment was accomplished as follows. Inspection of 
the data indicated that only 7% of the total loading is 
viscoelastic. This implied that 93% of the correction for the 
pi a tf orm drift is due to the change in the elastic component of 
the load. A cubic polynomial was least-squares fit to the stair- 
step hysteretic loading and unloading data. The resulting cubic 
curve is shown in Figure 4. It passes through the thin quasi- 
static hysteresis loop, and is an approximation to the elastic 
load-displacement curve. The platform location at the start of 
the load-relaxation, 2.194 in, is a reference position from which 
the drift can be measured. The total relaxing load was adjusted 
by the computed difference in the elastic component of the load 
due to the platform drifting. The adjusted load-relaxation curve 
is shown in Figure 5 . Note that the adjusted curve has more noise 
than the original curve. The noise is due to the fact that the 
pl a tform displacement data does not contain as many significant 
digits as the platform load data. However, the least— squares fit 
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to the noisy adjusted data presented in Figure 5 indicates that 
the resulting adjustment is significant. 

3 

With j = {10, 100, 1000}, the constrained least-squares 

Prony series for the adjusted load-relaxation data is 


-t 


-t 


m= 


1 + 0.01836 e }0 +0.01630 e 100 +0.03650 e 1000 


* 18523 


(15) 


More accurate least-squares curve fits can be obtained with a 
larger number of time constants (spaced more closely, etc.) 
However, the use of a large number of time constants slows down 
the finite element algorithm. All of the calculations below were 
performed employing the Prony series represented by Equation (15). 


Selection of Time Step 

Prior to running the viscous version of the finite element 
code a one-dimensional numerical simulation of the stair-step test 
was made. The one-dimensional equations are solved quickly and 
produce insight on the size of the time steps required in the 
finite element simulation. A schematic of the one-dimensional 
material model is shown in Figure 6. The quasi-static equilibrium 
equations for the model are 


3 

f(t) = fe(x) + k v'L a j X j 


n = 1 



dt 



dx 

dt 


for j= 1,2,3 


(16) 


where Tjj ^ are the Prony series coefficients, f e (x ) is the 

polynomial representation of the nonlinear elastic load 

(ofj + a 2 + cc 3 ) * 18523 lbs 

displacement curve, k v = is the total 

2.194 in 
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viscous stiffness obtained from the step strain relaxation data, 
and are the internal variables used to compute the viscous 


force . 


A plot of the measured stair-step displacement data is shown 
in Figure 7 . A piecewise linear representation of this curve was 
obtained by using straight lines between the corner points of the 
stair-steps. As noted above, close inspection of the data reveals 
that the platform was drifting during the relaxation intervals. 

The simulation of the stair-step loading that results from 
integrating Equations (16) is shown in Figure 8. Load- 
displacement plots at the peak load, computed with time steps of 
0 . 2 sec and 1 . 0 sec , are shown in Figure 9 . The plots for each of 
these two time steps agree well. Smaller time steps did not 
provide any additional information. The values of the time 

constants, {^„} n=] / in the differential equations indicate that a 


time step of 2.0 sec or larger (to integrate e 10 with the 
trapezoidal method) can be used. However, when a time step of 2.0 
sec or larger was used, the errors introduced by missing the 
details of the stair-step ramping action were too large to accept. 
A time step of 1.0 sec was selected to for the finite element 
computations . 


Finite Element Simulation 

The tire ' s finite element mesh and a sketch of the loading 
platform are shown in Figure 10. The mesh is similar to the 

*7 

"Model 1" mesh employed by Tanner. The elastic model has 540 

elements and 28,565 degrees of freedom (not including the Lagrange 
multipliers used for points that come into contact.) An 
additional 103,680 internal variables were added to program the 
solution algorithm for the material model described above. The 
platform surface is frictionless. Computed elastic and 

viscoelastic load-displacement curves, obtained by enforcing the 
stair-step tire rim displacement are shown in Figure 11. 
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A comparison of the curves in Figures 8 and 11 indicates that 
the finite element model is stiffer than the one-dimensional 
model. This is because the elastic component of the finite 
element model is stiffer than the elastic component of the one- 
dimensional model. Since the one-dimensional model closely 
represents the measured data, the finite element model employed 
here produces a load-displacement curve which is above the 
measured data curve. The finite element load-displacement 
hysteresis loop and the measured hysteresis loop are shown in 
Figure 12. The measured loop encloses more area than the computed 
loop. This indicates that the simulation underestimated the 
viscous energy lost during the test. 


CONCLUDING REMARKS 

An algorithm for converting elastic structural elements based 
on the mixed Hellinger-Resiner mixed variational principle to 
viscoelastic structural elements was presented. The thirteen node 

large displacement thick-shell element derived by Noor and 
4 

Hartley was employed to describe the algorithm. A finite element 

tire model based on this shell element, and used by Tanner ' ' to 

analyze tire footprints was modified so that the tire material 
would represent a Maxwell solid. Load-displacement data from a 
stair-step loading test was computationally simulated. The 
computed stair-step hysteresis loop indicated less viscous loss 
than the measured loop. The new computational algorithm 
functioned successfully. This algorithm can be applied to all 
structural elements of either displacement or mixed type. 
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Figure 1 . Shell displacement and force 
variables . 
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Figure 2. Shell finite element displacement 
and stress nodal variables. 
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Load on Tire (lbs) 


Load on Tire 



Time (sec) 


Figure 3 . Tire load-relaxation and 
rim displacement data. 
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Tire Rim Displacement (in) 



Load (lbs) 



Figure 4. Tire stair-step loading data. 
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Tire Load (lbs) 



Figure 5 . Adjusted tire load-relaxation 
data. 
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Figure 6. One dimensional material 

model. 
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Figure 8. Viscous and elastic loads for 
one-dimensional stair-step 
simulation. 
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Figure 9 . Computed tire loading for 
time steps 0.2 and 1.0 sec., 
one-dimensional simulation. 
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Rigid platform 



rim displacement 

Figure 10. Tire finite element model. 
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Figure 12. Finite element and test 

stair-step hysteresis curves . 
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